function [varargout]=Equ_HH(PP,SS,MODEL,JacobianFlag,Dev_Input,Dev_Output,TimePrintFlag)

%% Preliminaries
VAR         =   MODEL.VAR;
HH          =   MODEL.EquBlock.HH;

if nargin<=3
    JacobianFlag    =   0;
end
if nargin<=4
    Dev_Input       =   zeros(HH.Dim.In,1);
    Dev_Output      =   zeros(HH.Dim.Out,1);
elseif ( isempty(Dev_Input) || isempty(Dev_Output) )
    Dev_Input       =   zeros(HH.Dim.In,1);
    Dev_Output      =   zeros(HH.Dim.Out,1);
end
if nargin<=6
    TimePrintFlag   =   0;
end
if TimePrintFlag
    TimeStart       =   tic;
end
%% Evaluation
%--------------------------------------------------------------------------
% Pre-Transformation: Input Deviations to Input Levels
%--------------------------------------------------------------------------

for ii=1:HH.Num.In
    TempVar     =   HH.Var.In{ii};
    TempDev     =   Dev_Input(HH.LocIdx.In.(TempVar));
    TempLevel   =   VAR.(TempVar).Dev2Level(TempDev);
    eval(['Dev_',TempVar,'=TempDev;'])
    eval([TempVar,'=TempLevel;'])
end
Dev_Decision=   [Dev_Pir;Dev_T;Dev_Profit;Dev_TAU;Dev_ir;Dev_w;Dev_ValFunp];
Dev_Distribution ...
            =   [Dev_Dist];

%--------------------------------------------------------------------------
% 1. Decision Problem
%--------------------------------------------------------------------------
if max(abs(Dev_Decision))<eps && JacobianFlag>0
    ValFun      =   VAR.ValFun.SS;
    Policy      =   SS.DistPolicy;
    FlowVec     =   SS.TrProb.FlowVec;
    CoefMat     =   SS.TrProb.CoefMat;
    UnitTrProb  =   SS.TrProb.UnitTrMat;
else
    %----------------------------------------------------------------------
    % Solve the Decision Problem Sequentially
    %----------------------------------------------------------------------
    % Start Step: inverse State Reduction, ValFunp to EVCoefp
    UCoefp      =   SS.UCoef+ValFunp;
    % Step 1: from UCoefp to VbarCoef
    Opt_Vbar    =   Fun_Vbar(PP,SS,T,Profit,TAU,ir,w,UCoefp,SS.FunApp.V.State);
    VbarCoef_Mat=   SS.FunApp.V.BasMat_0\reshape(Opt_Vbar,SS.FunApp.V.MatSize_State);
    VbarCoef    =   VbarCoef_Mat(:);
    % Step 2: from VCoef to VbarCoef
    Opt_U       =   Fun_U(PP,SS,Pir,VbarCoef,SS.FunApp.U.State);
    UCoef_Mat   =   SS.FunApp.U.BasMat_0\reshape(Opt_U,SS.FunApp.U.MatSize_State);
    UCoef       =   UCoef_Mat(:);
    % End Step: UCoef to ValFun
    ValFun      =   UCoef-SS.UCoef;
    %----------------------------------------------------------------------
    % Solve the Policy over Distribution Grid and Transition Prob.
    %----------------------------------------------------------------------
    % Policy over the Histogram Grid
    Policy      =   Policy_Solver(PP,SS,T,Profit,TAU,ir,w,UCoefp,SS.DistApp.fz.State);
    % Transition Prob. Matrix determined by Policy
    [FlowVec,CoefMat,UnitTrProb] ...
                =   Distribution_TrProbMat(PP,SS,Pir,Policy);  
    
end

%--------------------------------------------------------------------------
% 2. Distribution Evolution
%--------------------------------------------------------------------------
if max(abs(Dev_Distribution))<eps && JacobianFlag>0
    QW_U        =   SS.Dist_fz.QW_fz;
else
    % inverse State Reduction: (Marginal_f,Marginal_zo) to (f,zo)
    QW_U        =   DimReduction_Distribution(PP,SS,'UnZip',Dist,SS.CopulaF);
end

% Evolution: U2Up
QW_Up       =   CoefMat*QW_U+FlowVec;
% State Reduction:
[Distp,~]   =   DimReduction_Distribution(PP,SS,'Zip',QW_Up);

%--------------------------------------------------------------------------
% 3. Aggregation
%--------------------------------------------------------------------------
Agg         =   Aggregation(PP,SS,QW_U,Policy,UnitTrProb);

B           =   Agg.B;
C           =   Agg.C;
L           =   Agg.L;

AggStat     =   [Agg.Avg.c.poor;Agg.Avg.c.rich];
%--------------------------------------------------------------------------
% Post-Transformation: Output Levels to Ouput Deviations
%--------------------------------------------------------------------------
Dev_Left    =   zeros(HH.Dim.Out,1);
for ii=1:HH.Num.Out
    TempVar     =   HH.Var.Out{ii};
    eval(['TempLevel=',TempVar,';']);
    TempDev     =   VAR.(TempVar).Level2Dev(TempLevel);
    Dev_Left(HH.LocIdx.Out.(TempVar)) ...
                =   TempDev;
end
%% Return Objectives

switch JacobianFlag
    case 0
        varargout{1}    =   -Dev_Output+Dev_Left;
    case 1
        TempFun         =   @(Dev_Input)Equ_HH(PP,SS,MODEL,0,Dev_Input,Dev_Output);
        Jac             =   struct();
        Jac.Input       =   parfdjac(TempFun,Dev_Input);
        Jac.Output      =   -speye(MODEL.EquBlock.HH.Dim.Out);
        
        for ii=1:HH.Num.In
            TempVar         =   HH.Var.In{ii};
            Jac.(TempVar)   =   Jac.Input(:,HH.LocIdx.In.(TempVar));
        end
        for ii=1:HH.Num.Out
            TempVar         =   HH.Var.Out{ii};
            Jac.(TempVar)   =   Jac.Output(:,HH.LocIdx.Out.(TempVar));
        end
        varargout{1}    =   Jac;
end
if TimePrintFlag
    TimeElapsed     =   toc(TimeStart);
    fprintf(['Jacobian of HH Block is done: ',...
             'Elapse Time=',num2str(TimeElapsed,'%.2f'),'s','\n']);
end

